{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "<center>\n",
    "<img src=\"https://raw.githubusercontent.com/Yorko/mlcourse.ai/master/img/ods_stickers.jpg\" />\n",
    "    \n",
    "## [mlcourse.ai](https://mlcourse.ai) - Open Machine Learning Course\n",
    "\n",
    "<center>\n",
    "Auteur: [Dmitriy Sergeyev](https://github.com/DmitrySerg), Data Scientist @Zeptolab, chargé de cours au Center of Mathematical Finance de MSU.  \n",
    "    Traduit par: @borowis et [Ousmane Cissé](https://fr.linkedin.com/in/ousmane-cisse).  \n",
    "Ce matériel est soumis aux termes et conditions de la licence [Creative Commons CC BY-NC-SA 4.0](https://creativecommons.org/licenses/by-nc-sa/4.0/).  \n",
    "L'utilisation gratuite est autorisée à des fins non commerciales."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "\n",
    "# <center> Analyse des séries temporelles en Python</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Bonjour à tous !\n",
    "\n",
    "Voyons comment travailler avec les séries chronologiques en Python: quelles méthodes et quels modèles pouvons-nous utiliser pour la prédiction, ce qu'est le lissage exponentiel double et triple, que faire si la stationnarité n'est pas votre chose préférée, comment construire SARIMA et rester en vie , comment faire des prédictions en utilisant xgboost ... De plus, tout cela sera appliqué à des exemples (difficiles) du monde réel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Plan de l'article:\n",
    "1. [Introduction](#Introduction)\n",
    "   - [Forecast quality metrics](#Forecast-quality-metrics)\n",
    "2. [Move, smoothe, evaluate](#Move,-smoothe,-evaluate)\n",
    "   - Rolling window estimations\n",
    "   - Exponential smoothing, Holt-Winters model\n",
    "   - Time-series cross validation, parameters selection\n",
    "3. [Econometric approach](#Econometric-approach)\n",
    "   - Stationarity, unit root\n",
    "   - Getting rid of non-stationarity\n",
    "   - SARIMA intuition and model building\n",
    "4. [Linear (and not quite) models for time series](#Linear-(and-not-quite)-models-for-time-series)\n",
    "   - [Feature extraction](#Feature-extraction)\n",
    "   - [Time series lags](#Time-series-lags)\n",
    "   - [Target encoding](#Target-encoding)\n",
    "   - [Regularization and feature selection](#Regularization-and-feature-selection)\n",
    "   - [Boosting](#Boosting)\n",
    "5. [Conclusion](#Conclusion)\n",
    "6. [Demo assignment](#Demo-assignment)\n",
    "7. [Useful resources](#Useful-resources)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Dans mon travail quotidien, je rencontre des tâches liées à des séries chronologiques presque tous les jours. Les questions les plus fréquemment posées sont les suivantes: que se passera-t-il avec nos mesures le lendemain / la semaine / le mois / etc., combien d'utilisateurs installeront notre application, combien de temps ils passeront en ligne, combien d'actions les utilisateurs effectueront-ils, etc. Nous pouvons aborder ces tâches de prédiction en utilisant différentes méthodes en fonction de la qualité requise de la prédiction, de la durée de la période de prévision et, bien sûr, du délai dans lequel nous devons choisir les caractéristiques et régler les paramètres pour obtenir les résultats souhaités."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Introduction\n",
    "\n",
    "Nous commençons par une simple [définition](https://en.wikipedia.org/wiki/Time_series) des séries chronologiques:\n",
    "> *Séries temporelles* est une série de points de données indexés (ou répertoriés ou représentés graphiquement) dans un ordre chronologique.\n",
    "\n",
    "Par conséquent, les données sont organisées par des horodatages relativement déterministes et peuvent, par rapport aux données d'échantillonnage aléatoire, contenir des informations supplémentaires que nous pouvons extraire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Importons quelques bibliothèques. Tout d'abord, nous aurons besoin de la bibliothèque [statsmodels](http://statsmodels.sourceforge.net/stable/), qui possède de nombreuses fonctions de modélisation statistique, y compris des séries chronologiques. Pour les afficionados de R qui ont dû passer à Python, `statsmodels` vous sera certainement plus familier car il prend en charge les définitions de modèles   \n",
    "comme «Salaire ~ Âge + Éducation»."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt  # plots\n",
    "import numpy as np  # vectors and matrices\n",
    "import pandas as pd  # tables and data manipulations\n",
    "import seaborn as sns  # more plots\n",
    "\n",
    "sns.set()\n",
    "\n",
    "import warnings  # `do not disturbe` mode\n",
    "from itertools import product  # some useful functions\n",
    "\n",
    "import scipy.stats as scs\n",
    "import statsmodels.api as sm\n",
    "import statsmodels.formula.api as smf  # statistics and econometrics\n",
    "import statsmodels.tsa.api as smt\n",
    "from dateutil.relativedelta import relativedelta  # working with dates with style\n",
    "from scipy.optimize import minimize  # for function minimization\n",
    "from tqdm import tqdm_notebook\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "À titre d'exemple, regardons les vraies données de jeux mobiles. Plus précisément, nous examinerons les annonces regardées par heure et les dépenses en devises dans le jeu par jour:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ads = pd.read_csv(\"../../data/ads.csv\", index_col=[\"Time\"], parse_dates=[\"Time\"])\n",
    "currency = pd.read_csv(\n",
    "    \"../../data/currency.csv\", index_col=[\"Time\"], parse_dates=[\"Time\"]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 6))\n",
    "plt.plot(ads.Ads)\n",
    "plt.title(\"Ads watched (hourly data)\")\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(15, 7))\n",
    "plt.plot(currency.GEMS_GEMS_SPENT)\n",
    "plt.title(\"In-game currency spent (daily data)\")\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Mesures de qualité des prévisions\n",
    "\n",
    "Avant de commencer les prévisions, comprenons comment mesurer la qualité de nos prévisions et jetons un œil aux mesures les plus couramment utilisées."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "- [R squared](http://scikit-learn.org/stable/modules/model_evaluation.html#r2-score-the-coefficient-of-determination): coefficient de détermination (en économétrie, cela peut être interprété comme le pourcentage de variance expliqué par le modèle), $(-\\infty, 1]$\n",
    "\n",
    "$R^2 = 1 - \\frac{SS_{res}}{SS_{tot}}$\n",
    "\n",
    "```python\n",
    "sklearn.metrics.r2_score\n",
    "```\n",
    "---\n",
    "- [Mean Absolute Error](http://scikit-learn.org/stable/modules/model_evaluation.html#mean-absolute-error): il s'agit d'une métrique interprétable car elle a la même unité de mesure que la série initiale, $[0, +\\infty)$\n",
    "\n",
    "$MAE = \\frac{\\sum\\limits_{i=1}^{n} |y_i - \\hat{y}_i|}{n}$\n",
    "\n",
    "```python\n",
    "sklearn.metrics.mean_absolute_error\n",
    "```\n",
    "---\n",
    "                        \n",
    "- [Median Absolute Error](http://scikit-learn.org/stable/modules/model_evaluation.html#median-absolute-error): encore une fois, une métrique interprétable qui est particulièrement intéressante car elle est robuste aux valeurs aberrantes, $[0, +\\infty)$\n",
    "\n",
    "$MedAE = median(|y_1 - \\hat{y}_1|, ... , |y_n - \\hat{y}_n|)$\n",
    "\n",
    "```python\n",
    "sklearn.metrics.median_absolute_error\n",
    "\n",
    "```\n",
    "---\n",
    "\n",
    "- [Mean Squared Error](http://scikit-learn.org/stable/modules/model_evaluation.html#mean-squared-error): la métrique la plus couramment utilisée qui donne une pénalité plus élevée aux grosses erreurs et vice versa, $[0, +\\infty)$\n",
    "\n",
    "$MSE = \\frac{1}{n}\\sum\\limits_{i=1}^{n} (y_i - \\hat{y}_i)^2$\n",
    "\n",
    "```python\n",
    "sklearn.metrics.mean_squared_error\n",
    "```\n",
    "---\n",
    "\n",
    "- [Mean Squared Logarithmic Error](http://scikit-learn.org/stable/modules/model_evaluation.html#mean-squared-logarithmic-error): pratiquement, la même chose que MSE, mais nous prenez le logarithme de la série. En conséquence, nous accordons également plus de poids aux petites erreurs. Ceci est généralement utilisé lorsque les données ont des tendances exponentielles, $[0, +\\infty)$\n",
    "\n",
    "$MSLE = \\frac{1}{n}\\sum\\limits_{i=1}^{n} (log(1+y_i) - log(1+\\hat{y}_i))^2$\n",
    "\n",
    "```python\n",
    "sklearn.metrics.mean_squared_log_error\n",
    "```\n",
    "---\n",
    "\n",
    "- Mean Absolute Percentage Error: c'est la même chose que MAE mais est calculé en pourcentage, ce qui est très pratique lorsque vous voulez expliquer la qualité du modèle à la direction, $[0, +\\infty)$\n",
    "\n",
    "$MAPE = \\frac{100}{n}\\sum\\limits_{i=1}^{n} \\frac{|y_i - \\hat{y}_i|}{y_i}$\n",
    "\n",
    "```python\n",
    "def mean_absolute_percentage_error(y_true, y_pred): \n",
    "    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100\n",
    "```\n",
    "</div><i class=\"fa fa-lightbulb-o \"></i>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importing everything from above\n",
    "\n",
    "from sklearn.metrics import (\n",
    "    mean_absolute_error,\n",
    "    mean_squared_error,\n",
    "    mean_squared_log_error,\n",
    "    median_absolute_error,\n",
    "    r2_score,\n",
    ")\n",
    "\n",
    "\n",
    "def mean_absolute_percentage_error(y_true, y_pred):\n",
    "    return np.mean(np.abs((y_true - y_pred) / y_true)) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant que nous savons comment mesurer la qualité des prévisions, voyons quelles mesures nous pouvons utiliser et comment traduire les résultats pour le Directeur. Après cela, un petit détail reste - la construction du modèle."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Déplacer, lisser, évaluer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Commençons par une hypothèse naïve: \"demain sera le même qu'aujourd'hui\". Cependant, au lieu d'un modèle comme $\\hat{y}_{t} = y_{t-1}$ (qui est en fait une excellente base de référence pour tout problème de prédiction de séries chronologiques et parfois impossible à battre), nous supposerons que la valeur future de notre variable dépend de la moyenne de ses valeurs précédentes $k$. Par conséquent, nous utiliserons la **moyenne mobile**.\n",
    "\n",
    "$\\hat{y}_{t} = \\frac{1}{k} \\displaystyle\\sum^{k}_{n=1} y_{t-n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def moving_average(series, n):\n",
    "    \"\"\"\n",
    "        Calculate average of last n observations\n",
    "    \"\"\"\n",
    "    return np.average(series[-n:])\n",
    "\n",
    "\n",
    "moving_average(ads, 24)  # prediction for the last observed day (past 24 hours)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Malheureusement, nous ne pouvons pas faire de prédictions dans le futur - pour obtenir la valeur de l'étape suivante, nous avons besoin que les valeurs précédentes soient réellement observées. Mais la moyenne mobile a un autre cas d'utilisation - lisser la série chronologique d'origine pour identifier les tendances. Pandas a une implémentation disponible avec [`DataFrame.rolling (window).mean ()`](http://pandas.pydata.org/pandas-docs/stable/generated/pandas.DataFrame.rolling.html). Plus la fenêtre est large, plus la tendance est fluide. Dans le cas de données très bruyantes, que l'on rencontre souvent en finance, cette procédure peut aider à détecter des schémas (patterns) courants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def plotMovingAverage(\n",
    "    series, window, plot_intervals=False, scale=1.96, plot_anomalies=False\n",
    "):\n",
    "\n",
    "    \"\"\"\n",
    "        series - dataframe with timeseries\n",
    "        window - rolling window size \n",
    "        plot_intervals - show confidence intervals\n",
    "        plot_anomalies - show anomalies \n",
    "\n",
    "    \"\"\"\n",
    "    rolling_mean = series.rolling(window=window).mean()\n",
    "\n",
    "    plt.figure(figsize=(15, 5))\n",
    "    plt.title(\"Moving average\\n window size = {}\".format(window))\n",
    "    plt.plot(rolling_mean, \"g\", label=\"Rolling mean trend\")\n",
    "\n",
    "    # Plot confidence intervals for smoothed values\n",
    "    if plot_intervals:\n",
    "        mae = mean_absolute_error(series[window:], rolling_mean[window:])\n",
    "        deviation = np.std(series[window:] - rolling_mean[window:])\n",
    "        lower_bond = rolling_mean - (mae + scale * deviation)\n",
    "        upper_bond = rolling_mean + (mae + scale * deviation)\n",
    "        plt.plot(upper_bond, \"r--\", label=\"Upper Bond / Lower Bond\")\n",
    "        plt.plot(lower_bond, \"r--\")\n",
    "\n",
    "        # Having the intervals, find abnormal values\n",
    "        if plot_anomalies:\n",
    "            anomalies = pd.DataFrame(index=series.index, columns=series.columns)\n",
    "            anomalies[series < lower_bond] = series[series < lower_bond]\n",
    "            anomalies[series > upper_bond] = series[series > upper_bond]\n",
    "            plt.plot(anomalies, \"ro\", markersize=10)\n",
    "\n",
    "    plt.plot(series[window:], label=\"Actual values\")\n",
    "    plt.legend(loc=\"upper left\")\n",
    "    plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Lissons les 4 dernières heures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(ads, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Essayons maintenant le lissage des 12 dernières heures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(ads, 12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant en lissant les dernières 24 heures, nous obtenons la tendance quotidienne."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(ads, 24)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Lorsque nous avons appliqué un lissage quotidien sur les données horaires, nous avons pu voir clairement la dynamique des annonces vues. Pendant les week-ends, les valeurs sont plus élevées (plus de temps pour jouer le week-end) tandis que moins d'annonces sont regardées en semaine."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous pouvons également tracer des intervalles de confiance pour nos valeurs lissées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(ads, 4, plot_intervals=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant, créons un système de détection d'anomalies simple à l'aide de la moyenne mobile. Malheureusement, dans cet ensemble de données particulier, tout est plus ou moins normal, nous allons donc intentionnellement créer des valeurs anormales dans notre dataframe `ads_anomaly`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ads_anomaly = ads.copy()\n",
    "ads_anomaly.iloc[-20] = ads_anomaly.iloc[-20] * 0.2  # say we have 80% drop of ads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voyons si cette méthode simple peut détecter l'anomalie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(ads_anomaly, 4, plot_intervals=True, plot_anomalies=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "Neat! What about the second series?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Et la deuxième série?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotMovingAverage(\n",
    "    currency, 7, plot_intervals=True, plot_anomalies=True\n",
    ")  # weekly smoothing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Oh non, ce n'était pas aussi génial! Ici, nous pouvons voir l'inconvénient de notre approche simple - elle n'a pas saisi la saisonnalité mensuelle dans nos données et a marqué presque tous les pics de 30 jours comme des anomalies. Si vous voulez éviter les faux positifs, il est préférable d'envisager des modèles plus complexes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "**Moyenne pondérée** (Weighted average) est une simple modification de la moyenne mobile. Les poids totalisent `1` avec des poids plus importants attribués à des observations plus récentes.\n",
    "\n",
    "\n",
    "$\\hat{y}_{t} = \\displaystyle\\sum^{k}_{n=1} \\omega_n y_{t+1-n}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def weighted_average(series, weights):\n",
    "    \"\"\"\n",
    "        Calculate weighted average on the series.\n",
    "        Assuming weights are sorted in descending order\n",
    "        (larger weights are assigned to more recent observations).\n",
    "    \"\"\"\n",
    "    result = 0.0\n",
    "    for n in range(len(weights)):\n",
    "        result += series.iloc[-n - 1] * weights[n]\n",
    "    return float(result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "weighted_average(ads, [0.6, 0.3, 0.1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# just checking\n",
    "0.6 * ads.iloc[-1] + 0.3 * ads.iloc[-2] + 0.1 * ads.iloc[-3]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Lissage exponentiel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voyons maintenant ce qui se passe si, au lieu de pondérer les dernières valeurs $k$ de la série chronologique, nous commençons à pondérer toutes les observations disponibles tout en diminuant de façon exponentielle les poids à mesure que nous remontons dans le temps. Il existe une formule pour **[lissage exponentiel](https://en.wikipedia.org/wiki/Exponential_smoothing)** qui nous aidera avec ceci:\n",
    "\n",
    "$$\\hat{y}_{t} = \\alpha \\cdot y_t + (1-\\alpha) \\cdot \\hat y_{t-1} $$\n",
    "\n",
    "Ici, la valeur du modèle est une moyenne pondérée entre la valeur réelle actuelle et les valeurs précédentes du modèle. Le poids $\\alpha$ est appelé facteur de lissage. Il définit la vitesse à laquelle nous \"oublierons\" la dernière observation vraie disponible. Plus $\\alpha$ est petit, plus les observations précédentes ont d'influence et plus la série est fluide.\n",
    "\n",
    "L'exponentialité est cachée dans la récursivité de la fonction - nous multiplions par $(1-\\alpha)$ à chaque fois, qui contient déjà une multiplication par $(1-\\alpha)$ des valeurs de modèle précédentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def exponential_smoothing(series, alpha):\n",
    "    \"\"\"\n",
    "        series - dataset with timestamps\n",
    "        alpha - float [0.0, 1.0], smoothing parameter\n",
    "    \"\"\"\n",
    "    result = [series[0]]  # first value is same as series\n",
    "    for n in range(1, len(series)):\n",
    "        result.append(alpha * series[n] + (1 - alpha) * result[n - 1])\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def plotExponentialSmoothing(series, alphas):\n",
    "    \"\"\"\n",
    "        Plots exponential smoothing with different alphas\n",
    "        \n",
    "        series - dataset with timestamps\n",
    "        alphas - list of floats, smoothing parameters\n",
    "        \n",
    "    \"\"\"\n",
    "    with plt.style.context(\"seaborn-white\"):\n",
    "        plt.figure(figsize=(15, 7))\n",
    "        for alpha in alphas:\n",
    "            plt.plot(\n",
    "                exponential_smoothing(series, alpha), label=\"Alpha {}\".format(alpha)\n",
    "            )\n",
    "        plt.plot(series.values, \"c\", label=\"Actual\")\n",
    "        plt.legend(loc=\"best\")\n",
    "        plt.axis(\"tight\")\n",
    "        plt.title(\"Exponential Smoothing\")\n",
    "        plt.grid(True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotExponentialSmoothing(ads.Ads, [0.3, 0.05])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotExponentialSmoothing(currency.GEMS_GEMS_SPENT, [0.3, 0.05])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Double lissage exponentiel"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Jusqu'à présent, les méthodes dont nous avons discuté, ont porté pour une seule prédiction de points futurs (avec un bon lissage). C'est cool, mais ce n'est pas suffisant non plus. Étendons le lissage exponentiel afin de pouvoir prédire deux points futurs (bien sûr, nous inclurons également plus de lissage).\n",
    "\n",
    "La décomposition en séries nous aidera - nous obtenons deux composantes: $\\ell$ d'interception (c'est-à-dire de niveau) et $b$ de pente (c'est-à-dire de tendance). Nous avons appris à prédire l'interception (ou la valeur de série attendue) avec nos méthodes précédentes; maintenant, nous appliquerons le même lissage exponentiel à la tendance en supposant que la direction future des changements de séries chronologiques dépend des changements pondérés précédents. En conséquence, nous obtenons l'ensemble de fonctions suivant:\n",
    "\n",
    "$$\\ell_x = \\alpha y_x + (1-\\alpha)(\\ell_{x-1} + b_{x-1})$$\n",
    "\n",
    "$$b_x = \\beta(\\ell_x - \\ell_{x-1}) + (1-\\beta)b_{x-1}$$\n",
    "\n",
    "$$\\hat{y}_{x+1} = \\ell_x + b_x$$\n",
    "\n",
    "La première décrit l'ordonnée à l'origine, qui, comme précédemment, dépend de la valeur actuelle de la série. Le deuxième terme est maintenant divisé en valeurs précédentes du niveau et de la tendance. La deuxième fonction décrit la tendance, qui dépend des changements de niveau à l'étape actuelle et de la valeur précédente de la tendance. Dans ce cas, le coefficient $\\beta$ est un poids pour le lissage exponentiel. La prédiction finale est la somme des valeurs du modèle de l'interception et de la tendance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     20
    ]
   },
   "outputs": [],
   "source": [
    "def double_exponential_smoothing(series, alpha, beta):\n",
    "    \"\"\"\n",
    "        series - dataset with timeseries\n",
    "        alpha - float [0.0, 1.0], smoothing parameter for level\n",
    "        beta - float [0.0, 1.0], smoothing parameter for trend\n",
    "    \"\"\"\n",
    "    # first value is same as series\n",
    "    result = [series[0]]\n",
    "    for n in range(1, len(series) + 1):\n",
    "        if n == 1:\n",
    "            level, trend = series[0], series[1] - series[0]\n",
    "        if n >= len(series):  # forecasting\n",
    "            value = result[-1]\n",
    "        else:\n",
    "            value = series[n]\n",
    "        last_level, level = level, alpha * value + (1 - alpha) * (level + trend)\n",
    "        trend = beta * (level - last_level) + (1 - beta) * trend\n",
    "        result.append(level + trend)\n",
    "    return result\n",
    "\n",
    "\n",
    "def plotDoubleExponentialSmoothing(series, alphas, betas):\n",
    "    \"\"\"\n",
    "        Plots double exponential smoothing with different alphas and betas\n",
    "        \n",
    "        series - dataset with timestamps\n",
    "        alphas - list of floats, smoothing parameters for level\n",
    "        betas - list of floats, smoothing parameters for trend\n",
    "    \"\"\"\n",
    "\n",
    "    with plt.style.context(\"seaborn-white\"):\n",
    "        plt.figure(figsize=(20, 8))\n",
    "        for alpha in alphas:\n",
    "            for beta in betas:\n",
    "                plt.plot(\n",
    "                    double_exponential_smoothing(series, alpha, beta),\n",
    "                    label=\"Alpha {}, beta {}\".format(alpha, beta),\n",
    "                )\n",
    "        plt.plot(series.values, label=\"Actual\")\n",
    "        plt.legend(loc=\"best\")\n",
    "        plt.axis(\"tight\")\n",
    "        plt.title(\"Double Exponential Smoothing\")\n",
    "        plt.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotDoubleExponentialSmoothing(ads.Ads, alphas=[0.9, 0.02], betas=[0.9, 0.02])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotDoubleExponentialSmoothing(\n",
    "    currency.GEMS_GEMS_SPENT, alphas=[0.9, 0.02], betas=[0.9, 0.02]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant, nous devons régler deux paramètres: $\\alpha$ et $\\beta$. Le premier est responsable du lissage de la série autour de la tendance, le second du lissage de la tendance elle-même. Plus les valeurs sont élevées, plus les observations les plus récentes auront de poids et moins la série de modèles sera lissée. Certaines combinaisons de paramètres peuvent produire des résultats étranges, surtout si elles sont définies manuellement. Nous nous pencherons dans un instant sur le choix automatique des paramètres ; avant cela, parlons du lissage triple exponentiel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Triple lissage exponentiel alias Holt-Winters\n",
    "\n",
    "Nous avons examiné le lissage exponentiel et le double lissage exponentiel. Cette fois, nous allons passer au _triple_ lissage exponentiel.\n",
    "\n",
    "Comme vous auriez pu le deviner, l'idée est d'ajouter une troisième composante - la saisonnalité. Cela signifie que nous ne devrions pas utiliser cette méthode si notre série chronologique ne devrait pas avoir de saisonnalité. Les composantes saisonnières du modèle expliqueront les variations répétées autour de l'interception et de la tendance, et elles seront spécifiées par la durée de la saison, en d'autres termes par la période après laquelle les variations se répètent. Pour chaque observation de la saison, il y a une composante distincte; par exemple, si la durée de la saison est de 7 jours (une saisonnalité hebdomadaire), nous aurons 7 composantes saisonnières, une pour chaque jour de la semaine.\n",
    "\n",
    "Avec cela, écrivons un nouveau système d'équations:\n",
    "\n",
    "$$\\ell_x = \\alpha(y_x - s_{x-L}) + (1-\\alpha)(\\ell_{x-1} + b_{x-1})$$\n",
    "\n",
    "$$b_x = \\beta(\\ell_x - \\ell_{x-1}) + (1-\\beta)b_{x-1}$$\n",
    "\n",
    "$$s_x = \\gamma(y_x - \\ell_x) + (1-\\gamma)s_{x-L}$$\n",
    "\n",
    "$$\\hat{y}_{x+m} = \\ell_x + mb_x + s_{x-L+1+(m-1)modL}$$\n",
    "\n",
    "L'interception dépend maintenant de la valeur actuelle de la série moins toute composante saisonnière correspondante. La tendance reste inchangée et la composante saisonnière dépend de la valeur actuelle de la série moins l'ordonnée à l'origine et de la valeur précédente de la composante. Prendre en compte que la composante est lissée à travers toutes les saisons disponibles; par exemple, si nous avons une composante lundi, elle ne sera calculée qu'en moyenne avec les autres lundis. Vous pouvez en savoir plus sur le fonctionnement de la moyenne et comment se fait l'approximation initiale de la tendance et des composantes saisonnières [ici](http://www.itl.nist.gov/div898/handbook/pmc/section4/pmc435.htm). Maintenant que nous avons la composante saisonnière, nous pouvons prédire non seulement un ou deux pas en avant, mais un futur arbitraire $m$, ce qui est très encourageant.\n",
    "\n",
    "Vous trouverez ci-dessous le code d'un modèle de lissage exponentiel triple, également connu sous les noms de famille de ses créateurs, Charles Holt et son élève Peter Winters. De plus, la méthode Brutlag a été incluse dans le modèle pour produire des intervalles de confiance:\n",
    "\n",
    "$$\\hat y_{max_x}=\\ell_{x−1}+b_{x−1}+s_{x−T}+m⋅d_{t−T}$$\n",
    "\n",
    "$$\\hat y_{min_x}=\\ell_{x−1}+b_{x−1}+s_{x−T}-m⋅d_{t−T}$$\n",
    "\n",
    "$$d_t=\\gamma∣y_t−\\hat y_t∣+(1−\\gamma)d_{t−T},$$\n",
    "\n",
    "où $T$ est la durée de la saison, $d$ est l'écart prévu. D'autres paramètres ont été tirés du triple lissage exponentiel. Vous pouvez en savoir plus sur la méthode et son applicabilité à la détection d'anomalies dans des séries chronologiques [ici](http://fedcsis.org/proceedings/2012/pliks/118.pdf)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "class HoltWinters:\n",
    "\n",
    "    \"\"\"\n",
    "    Holt-Winters model with the anomalies detection using Brutlag method\n",
    "    \n",
    "    # series - initial time series\n",
    "    # slen - length of a season\n",
    "    # alpha, beta, gamma - Holt-Winters model coefficients\n",
    "    # n_preds - predictions horizon\n",
    "    # scaling_factor - sets the width of the confidence interval by Brutlag (usually takes values from 2 to 3)\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, series, slen, alpha, beta, gamma, n_preds, scaling_factor=1.96):\n",
    "        self.series = series\n",
    "        self.slen = slen\n",
    "        self.alpha = alpha\n",
    "        self.beta = beta\n",
    "        self.gamma = gamma\n",
    "        self.n_preds = n_preds\n",
    "        self.scaling_factor = scaling_factor\n",
    "\n",
    "    def initial_trend(self):\n",
    "        sum = 0.0\n",
    "        for i in range(self.slen):\n",
    "            sum += float(self.series[i + self.slen] - self.series[i]) / self.slen\n",
    "        return sum / self.slen\n",
    "\n",
    "    def initial_seasonal_components(self):\n",
    "        seasonals = {}\n",
    "        season_averages = []\n",
    "        n_seasons = int(len(self.series) / self.slen)\n",
    "        # let's calculate season averages\n",
    "        for j in range(n_seasons):\n",
    "            season_averages.append(\n",
    "                sum(self.series[self.slen * j : self.slen * j + self.slen])\n",
    "                / float(self.slen)\n",
    "            )\n",
    "        # let's calculate initial values\n",
    "        for i in range(self.slen):\n",
    "            sum_of_vals_over_avg = 0.0\n",
    "            for j in range(n_seasons):\n",
    "                sum_of_vals_over_avg += (\n",
    "                    self.series[self.slen * j + i] - season_averages[j]\n",
    "                )\n",
    "            seasonals[i] = sum_of_vals_over_avg / n_seasons\n",
    "        return seasonals\n",
    "\n",
    "    def triple_exponential_smoothing(self):\n",
    "        self.result = []\n",
    "        self.Smooth = []\n",
    "        self.Season = []\n",
    "        self.Trend = []\n",
    "        self.PredictedDeviation = []\n",
    "        self.UpperBond = []\n",
    "        self.LowerBond = []\n",
    "\n",
    "        seasonals = self.initial_seasonal_components()\n",
    "\n",
    "        for i in range(len(self.series) + self.n_preds):\n",
    "            if i == 0:  # components initialization\n",
    "                smooth = self.series[0]\n",
    "                trend = self.initial_trend()\n",
    "                self.result.append(self.series[0])\n",
    "                self.Smooth.append(smooth)\n",
    "                self.Trend.append(trend)\n",
    "                self.Season.append(seasonals[i % self.slen])\n",
    "\n",
    "                self.PredictedDeviation.append(0)\n",
    "\n",
    "                self.UpperBond.append(\n",
    "                    self.result[0] + self.scaling_factor * self.PredictedDeviation[0]\n",
    "                )\n",
    "\n",
    "                self.LowerBond.append(\n",
    "                    self.result[0] - self.scaling_factor * self.PredictedDeviation[0]\n",
    "                )\n",
    "                continue\n",
    "\n",
    "            if i >= len(self.series):  # predicting\n",
    "                m = i - len(self.series) + 1\n",
    "                self.result.append((smooth + m * trend) + seasonals[i % self.slen])\n",
    "\n",
    "                # when predicting we increase uncertainty on each step\n",
    "                self.PredictedDeviation.append(self.PredictedDeviation[-1] * 1.01)\n",
    "\n",
    "            else:\n",
    "                val = self.series[i]\n",
    "                last_smooth, smooth = (\n",
    "                    smooth,\n",
    "                    self.alpha * (val - seasonals[i % self.slen])\n",
    "                    + (1 - self.alpha) * (smooth + trend),\n",
    "                )\n",
    "                trend = self.beta * (smooth - last_smooth) + (1 - self.beta) * trend\n",
    "                seasonals[i % self.slen] = (\n",
    "                    self.gamma * (val - smooth)\n",
    "                    + (1 - self.gamma) * seasonals[i % self.slen]\n",
    "                )\n",
    "                self.result.append(smooth + trend + seasonals[i % self.slen])\n",
    "\n",
    "                # Deviation is calculated according to Brutlag algorithm.\n",
    "                self.PredictedDeviation.append(\n",
    "                    self.gamma * np.abs(self.series[i] - self.result[i])\n",
    "                    + (1 - self.gamma) * self.PredictedDeviation[-1]\n",
    "                )\n",
    "\n",
    "            self.UpperBond.append(\n",
    "                self.result[-1] + self.scaling_factor * self.PredictedDeviation[-1]\n",
    "            )\n",
    "\n",
    "            self.LowerBond.append(\n",
    "                self.result[-1] - self.scaling_factor * self.PredictedDeviation[-1]\n",
    "            )\n",
    "\n",
    "            self.Smooth.append(smooth)\n",
    "            self.Trend.append(trend)\n",
    "            self.Season.append(seasonals[i % self.slen])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Validation croisée des séries chronologiques\n",
    "\n",
    "Avant de commencer à construire un modèle, voyons d'abord comment estimer automatiquement les paramètres du modèle.\n",
    "\n",
    "Il n'y a rien d'inhabituel ici; comme toujours, nous devons choisir une fonction de perte adaptée à la tâche qui nous dira à quel point le modèle se rapproche des données. Ensuite, en utilisant la validation croisée, nous évaluerons notre fonction de perte choisie pour les paramètres de modèle donnés, calculerons le gradient, ajusterons les paramètres de modèle, etc., pour finalement descendre au minimum global.\n",
    "\n",
    "Vous vous demandez peut-être comment effectuer la validation croisée pour les séries temporelles car les séries temporelles ont cette structure temporelle et on ne peut pas mélanger aléatoirement les valeurs dans un pli tout en préservant cette structure. Avec la randomisation, toutes les dépendances temporelles entre les observations seront perdues. C'est pourquoi nous devrons utiliser une approche plus délicate pour optimiser les paramètres du modèle. Je ne sais pas s'il y a un nom officiel à cela, mais sur [CrossValidated](https://stats.stackexchange.com/questions/14099/using-k-fold-cross-validation-for-time-series-model -sélection), où l'on peut trouver toutes les réponses sauf la réponse à la question ultime de la vie, de l'univers et de tout, le nom proposé pour cette méthode est \"validation croisée sur une base continue\".\n",
    "\n",
    "L'idée est assez simple - nous formons notre modèle sur un petit segment de la série chronologique du début jusqu'à certains $t$, faisons des prédictions pour les prochaines étapes de $t+n$ et calculons une erreur. Ensuite, nous étendons notre échantillon d'entraînement à la valeur $t+n$, faisons des prédictions de $t+n$ à $t+2*n$ et continuons de déplacer notre segment de test de la série chronologique jusqu'à ce que nous atteignions la dernière observation disponible. Par conséquent, nous avons autant de plis que $n$ se situera entre l'échantillon d'apprentissage initial et la dernière observation.\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/Yorko/mlcourse.ai/master/img/time_series_cv.png\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant, sachant comment configurer la validation croisée, nous pouvons trouver les paramètres optimaux pour le modèle Holt-Winters. Rappelons que nous avons une saisonnalité quotidienne dans les publicités, d'où le paramètre `slen = 24`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     2
    ]
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import TimeSeriesSplit  # you have everything done for you\n",
    "\n",
    "\n",
    "def timeseriesCVscore(params, series, loss_function=mean_squared_error, slen=24):\n",
    "    \"\"\"\n",
    "        Returns error on CV  \n",
    "        \n",
    "        params - vector of parameters for optimization\n",
    "        series - dataset with timeseries\n",
    "        slen - season length for Holt-Winters model\n",
    "    \"\"\"\n",
    "    # errors array\n",
    "    errors = []\n",
    "\n",
    "    values = series.values\n",
    "    alpha, beta, gamma = params\n",
    "\n",
    "    # set the number of folds for cross-validation\n",
    "    tscv = TimeSeriesSplit(n_splits=3)\n",
    "\n",
    "    # iterating over folds, train model on each, forecast and calculate error\n",
    "    for train, test in tscv.split(values):\n",
    "\n",
    "        model = HoltWinters(\n",
    "            series=values[train],\n",
    "            slen=slen,\n",
    "            alpha=alpha,\n",
    "            beta=beta,\n",
    "            gamma=gamma,\n",
    "            n_preds=len(test),\n",
    "        )\n",
    "        model.triple_exponential_smoothing()\n",
    "\n",
    "        predictions = model.result[-len(test) :]\n",
    "        actual = values[test]\n",
    "        error = loss_function(predictions, actual)\n",
    "        errors.append(error)\n",
    "\n",
    "    return np.mean(np.array(errors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Dans le modèle Holt-Winters, ainsi que dans les autres modèles de lissage exponentiel, il existe une contrainte sur la taille des paramètres de lissage, chacun variant de 0 à 1. Par conséquent, afin de minimiser notre fonction de perte, nous devons choisir un algorithme qui prend en charge les contraintes sur les paramètres du modèle. Dans notre cas, nous utiliserons le gradient conjugué de Newton tronqué."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "data = ads.Ads[:-20]  # leave some data for testing\n",
    "\n",
    "# initializing model parameters alpha, beta and gamma\n",
    "x = [0, 0, 0]\n",
    "\n",
    "# Minimizing the loss function\n",
    "opt = minimize(\n",
    "    timeseriesCVscore,\n",
    "    x0=x,\n",
    "    args=(data, mean_squared_log_error),\n",
    "    method=\"TNC\",\n",
    "    bounds=((0, 1), (0, 1), (0, 1)),\n",
    ")\n",
    "\n",
    "# Take optimal values...\n",
    "alpha_final, beta_final, gamma_final = opt.x\n",
    "print(alpha_final, beta_final, gamma_final)\n",
    "\n",
    "# ...and train the model with them, forecasting for the next 50 hours\n",
    "model = HoltWinters(\n",
    "    data,\n",
    "    slen=24,\n",
    "    alpha=alpha_final,\n",
    "    beta=beta_final,\n",
    "    gamma=gamma_final,\n",
    "    n_preds=50,\n",
    "    scaling_factor=3,\n",
    ")\n",
    "model.triple_exponential_smoothing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Ajoutons du code pour afficher les tracés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def plotHoltWinters(series, plot_intervals=False, plot_anomalies=False):\n",
    "    \"\"\"\n",
    "        series - dataset with timeseries\n",
    "        plot_intervals - show confidence intervals\n",
    "        plot_anomalies - show anomalies \n",
    "    \"\"\"\n",
    "\n",
    "    plt.figure(figsize=(20, 10))\n",
    "    plt.plot(model.result, label=\"Model\")\n",
    "    plt.plot(series.values, label=\"Actual\")\n",
    "    error = mean_absolute_percentage_error(series.values, model.result[: len(series)])\n",
    "    plt.title(\"Mean Absolute Percentage Error: {0:.2f}%\".format(error))\n",
    "\n",
    "    if plot_anomalies:\n",
    "        anomalies = np.array([np.NaN] * len(series))\n",
    "        anomalies[series.values < model.LowerBond[: len(series)]] = series.values[\n",
    "            series.values < model.LowerBond[: len(series)]\n",
    "        ]\n",
    "        anomalies[series.values > model.UpperBond[: len(series)]] = series.values[\n",
    "            series.values > model.UpperBond[: len(series)]\n",
    "        ]\n",
    "        plt.plot(anomalies, \"o\", markersize=10, label=\"Anomalies\")\n",
    "\n",
    "    if plot_intervals:\n",
    "        plt.plot(model.UpperBond, \"r--\", alpha=0.5, label=\"Up/Low confidence\")\n",
    "        plt.plot(model.LowerBond, \"r--\", alpha=0.5)\n",
    "        plt.fill_between(\n",
    "            x=range(0, len(model.result)),\n",
    "            y1=model.UpperBond,\n",
    "            y2=model.LowerBond,\n",
    "            alpha=0.2,\n",
    "            color=\"grey\",\n",
    "        )\n",
    "\n",
    "    plt.vlines(\n",
    "        len(series),\n",
    "        ymin=min(model.LowerBond),\n",
    "        ymax=max(model.UpperBond),\n",
    "        linestyles=\"dashed\",\n",
    "    )\n",
    "    plt.axvspan(len(series) - 20, len(model.result), alpha=0.3, color=\"lightgrey\")\n",
    "    plt.grid(True)\n",
    "    plt.axis(\"tight\")\n",
    "    plt.legend(loc=\"best\", fontsize=13);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotHoltWinters(ads.Ads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotHoltWinters(ads.Ads, plot_intervals=True, plot_anomalies=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "À en juger par les graphiques, notre modèle a réussi à approximer avec succès la série chronologique initiale, capturant la saisonnalité quotidienne, la tendance globale à la baisse et même certaines anomalies. Si vous regardez les écarts du modèle, vous pouvez clairement voir que le modèle réagit assez fortement aux changements dans la structure de la série, puis ramène rapidement l'écart aux valeurs normales, \"oubliant\" essentiellement le passé. Cette caractéristique du modèle nous permet de construire rapidement des systèmes de détection d'anomalies, même pour les données de série bruyantes, sans dépenser trop de temps et d'argent pour préparer les données et entraîner le modèle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(25, 5))\n",
    "plt.plot(model.PredictedDeviation)\n",
    "plt.grid(True)\n",
    "plt.axis(\"tight\")\n",
    "plt.title(\"Brutlag's predicted deviation\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous appliquerons le même algorithme pour la deuxième série qui, comme vous vous en souvenez peut-être, a une tendance et une saisonnalité de 30 jours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "data = currency.GEMS_GEMS_SPENT[:-50]\n",
    "slen = 30  # 30-day seasonality\n",
    "\n",
    "x = [0, 0, 0]\n",
    "\n",
    "opt = minimize(\n",
    "    timeseriesCVscore,\n",
    "    x0=x,\n",
    "    args=(data, mean_absolute_percentage_error, slen),\n",
    "    method=\"TNC\",\n",
    "    bounds=((0, 1), (0, 1), (0, 1)),\n",
    ")\n",
    "\n",
    "alpha_final, beta_final, gamma_final = opt.x\n",
    "print(alpha_final, beta_final, gamma_final)\n",
    "\n",
    "model = HoltWinters(\n",
    "    data,\n",
    "    slen=slen,\n",
    "    alpha=alpha_final,\n",
    "    beta=beta_final,\n",
    "    gamma=gamma_final,\n",
    "    n_preds=100,\n",
    "    scaling_factor=3,\n",
    ")\n",
    "model.triple_exponential_smoothing()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotHoltWinters(currency.GEMS_GEMS_SPENT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Cela semble bon! Le modèle a pris à la fois la tendance à la hausse et les pointes saisonnières et correspond assez bien aux données."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotHoltWinters(currency.GEMS_GEMS_SPENT, plot_intervals=True, plot_anomalies=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(20, 5))\n",
    "plt.plot(model.PredictedDeviation)\n",
    "plt.grid(True)\n",
    "plt.axis(\"tight\")\n",
    "plt.title(\"Brutlag's predicted deviation\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Approche économétrique"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "### Stationnarité (Stationarity)\n",
    "\n",
    "Avant de commencer la modélisation, nous devons mentionner une propriété importante des séries chronologiques: [**stationarity**](https://fr.wikipedia.org/wiki/Processus_stationnaire).\n",
    "\n",
    "Si un processus est stationnaire, cela signifie qu'il ne change pas ses propriétés statistiques dans le temps, à savoir sa moyenne et sa variance. (La constance de la variance est appelée [homoscédasticité](https://fr.wikipedia.org/wiki/Homosc%C3%A9dasticit%C3%A9)) La fonction de covariance ne dépend pas du temps; cela ne devrait dépendre que de la distance entre les observations. Vous pouvez le voir visuellement sur les images dans le post de [Sean Abu](http://www.seanabu.com/2016/03/22/time-series-seasonal-ARIMA-model-in-python/):\n",
    "\n",
    "- Le graphique rouge ci-dessous n'est pas stationnaire car la moyenne augmente avec le temps.\n",
    "\n",
    "<img src=\"https://habrastorage.org/files/20c/9d8/a63/20c9d8a633ec436f91dccd4aedcc6940.png\"/>\n",
    "\n",
    "- Nous n'avons pas eu de chance avec la variance et nous voyons la dispersion variable des valeurs dans le temps\n",
    "\n",
    "<img src=\"https://habrastorage.org/files/b88/eec/a67/b88eeca676d642449cab135273fd5a95.png\"/>\n",
    "\n",
    "- Enfin, la covariance du i ème terme et du (i + m) ème terme ne doit pas être fonction du temps. Dans le graphique suivant, vous remarquerez que l'écart se rapproche à mesure que le temps augmente. Par conséquent, la covariance n'est pas constante avec le temps dans le graphique de droite.\n",
    "\n",
    "<img src=\"https://habrastorage.org/files/2f6/1ee/cb2/2f61eecb20714352840748b826e38680.png\"/>\n",
    "\n",
    "Alors pourquoi la stationnarité est-elle si importante? Parce qu'il est facile de faire des prédictions sur une série stationnaire car on peut supposer que les futures propriétés statistiques ne seront pas différentes de celles actuellement observées. La plupart des modèles de séries chronologiques, d'une manière ou d'une autre, tentent de prédire ces propriétés (moyenne ou variance, par exemple). De fausses prédictions seraient fausses si la série originale n'était pas stationnaire. Malheureusement, la plupart des séries chronologiques que nous voyons en dehors des manuels scolaires ne sont pas fixes, mais nous pouvons (et devons) changer cela.\n",
    "\n",
    "Donc, pour lutter contre la non-stationnarité, nous devons pour ainsi dire connaître notre ennemi. Voyons comment nous pouvons le détecter. Nous examinerons le bruit blanc et les promenades aléatoires pour apprendre à passer gratuitement de l'un à l'autre."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Graphique du bruit blanc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "white_noise = np.random.normal(size=1000)\n",
    "with plt.style.context(\"bmh\"):\n",
    "    plt.figure(figsize=(15, 5))\n",
    "    plt.plot(white_noise)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Le processus généré par la distribution normale standard est stationnaire et oscille autour de 0 avec un écart de 1. Maintenant, sur la base de ce processus, nous allons générer un nouveau où chaque valeur suivante dépendra de la précédente: $x_t = \\rho x_{t-1} + e_t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voici le code pour afficher les tracés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def plotProcess(n_samples=1000, rho=0):\n",
    "    x = w = np.random.normal(size=n_samples)\n",
    "    for t in range(n_samples):\n",
    "        x[t] = rho * x[t - 1] + w[t]\n",
    "\n",
    "    with plt.style.context(\"bmh\"):\n",
    "        plt.figure(figsize=(10, 3))\n",
    "        plt.plot(x)\n",
    "        plt.title(\n",
    "            \"Rho {}\\n Dickey-Fuller p-value: {}\".format(\n",
    "                rho, round(sm.tsa.stattools.adfuller(x)[1], 3)\n",
    "            )\n",
    "        )\n",
    "\n",
    "\n",
    "for rho in [0, 0.6, 0.9, 1]:\n",
    "    plotProcess(rho=rho)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Sur le premier tracé, vous pouvez voir le même bruit blanc stationnaire qu'auparavant. Sur le deuxième graphique avec $\\rho$ augmenté à 0,6, des cycles plus larges sont apparus, mais il semble toujours stationnaire dans l'ensemble. Le troisième graphique s'écarte encore plus de la moyenne de 0 mais oscille toujours autour de la moyenne. Enfin, avec $\\rho=1$, nous avons un processus de marche aléatoire, c'est-à-dire une série chronologique non stationnaire.\n",
    "\n",
    "Cela se produit car, après avoir atteint la valeur critique, la série $x_t = \\rho x_{t-1} + e_t$ ne revient pas à sa valeur moyenne. Si nous soustrayons $x_{t-1}$ des deux côtés, nous obtiendrons $x_t - x_{t-1} = (\\rho - 1) x_{t-1} + e_t$, où l'expression de gauche est appelée la première différence. Si $\\rho=1$, alors la première différence nous donne $e_t$ de bruit blanc stationnaire. C'est l'idée principale derrière le [test de Dickey-Fuller](https://en.wikipedia.org/wiki/Dickey%E2%80%93Fuller_test) pour la stationnarité des séries chronologiques (test de la présence d'une racine unitaire). Si nous pouvons obtenir une série stationnaire à partir d'une série non stationnaire en utilisant la première différence, nous appelons ces séries intégrées d'ordre 1. L'hypothèse nulle du test est que la série chronologique est non stationnaire, ce qui a été rejeté sur les trois premiers graphiques et finalement accepté sur le dernier. Nous devons dire que la première différence n'est pas toujours suffisante pour obtenir une série stationnaire car le processus peut être intégré d'ordre d, d> 1 (et avoir plusieurs racines unitaires). Dans ce cas, on utilise le test Dickey-Fuller amélioré, qui vérifie plusieurs retards à la fois\n",
    "\n",
    "Nous pouvons lutter contre la non-stationnarité en utilisant différentes approches: diverses différences d'ordre, suppression de tendance et de saisonnalité, lissage et transformations comme Box-Cox ou logarithmique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Se débarrasser de la non-stationnarité et créer SARIMA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Construisons un modèle ARIMA en parcourant tous les étapes de création d'une série stationnaire."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Voici le code pour afficher les tracés."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def tsplot(y, lags=None, figsize=(12, 7), style=\"bmh\"):\n",
    "    \"\"\"\n",
    "        Plot time series, its ACF and PACF, calculate Dickey–Fuller test\n",
    "        \n",
    "        y - timeseries\n",
    "        lags - how many lags to include in ACF, PACF calculation\n",
    "    \"\"\"\n",
    "    if not isinstance(y, pd.Series):\n",
    "        y = pd.Series(y)\n",
    "\n",
    "    with plt.style.context(style):\n",
    "        fig = plt.figure(figsize=figsize)\n",
    "        layout = (2, 2)\n",
    "        ts_ax = plt.subplot2grid(layout, (0, 0), colspan=2)\n",
    "        acf_ax = plt.subplot2grid(layout, (1, 0))\n",
    "        pacf_ax = plt.subplot2grid(layout, (1, 1))\n",
    "\n",
    "        y.plot(ax=ts_ax)\n",
    "        p_value = sm.tsa.stattools.adfuller(y)[1]\n",
    "        ts_ax.set_title(\n",
    "            \"Time Series Analysis Plots\\n Dickey-Fuller: p={0:.5f}\".format(p_value)\n",
    "        )\n",
    "        smt.graphics.plot_acf(y, lags=lags, ax=acf_ax)\n",
    "        smt.graphics.plot_pacf(y, lags=lags, ax=pacf_ax)\n",
    "        plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsplot(ads.Ads, lags=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "_this outlier on partial autocorrelation plot looks like a statsmodels bug, partial autocorrelation shall be <= 1 like any correlation._\n",
    "\n",
    "Surprisingly, the initial series are stationary; the Dickey-Fuller test rejected the null hypothesis that a unit root is present. Actually, we can see this on the plot itself – we do not have a visible trend, so the mean is constant and the variance is pretty much stable. The only thing left is seasonality, which we have to deal with prior to modeling. To do so, let's take the \"seasonal difference\", which means a simple subtraction of the series from itself with a lag that equals the seasonal period."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "_cette valeur aberrante sur le tracé d'autocorrélation partielle ressemble à un bug des modèles statistiques, l'autocorrélation partielle doit être <= 1 comme toute corrélation._\n",
    "\n",
    "Étonnamment, les premières séries sont stationnaires; le test de Dickey-Fuller a rejeté l'hypothèse nulle selon laquelle une racine unitaire est présente. En fait, nous pouvons le voir sur le graphique lui-même - nous n'avons pas de tendance visible, donc la moyenne est constante et la variance est à peu près stable. La seule chose qui reste est la saisonnalité, à laquelle nous devons faire face avant la modélisation. Pour ce faire, prenons la «différence saisonnière», ce qui signifie une simple soustraction de la série d'elle-même avec un décalage égal à la période saisonnière."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ads_diff = ads.Ads - ads.Ads.shift(24)\n",
    "tsplot(ads_diff[24:], lags=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "C'est maintenant beaucoup mieux avec la saisonnalité visible disparue. Cependant, la fonction d'autocorrélation a encore trop de retards importants. Pour les supprimer, nous prendrons les premières différences, en soustrayant la série d'elle-même avec le décalage 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ads_diff = ads_diff - ads_diff.shift(1)\n",
    "tsplot(ads_diff[24 + 1 :], lags=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Parfait! Notre série ressemble maintenant à quelque chose d'indescriptible, oscillant autour de zéro. Le test de Dickey-Fuller indique qu'il est stationnaire et le nombre de pics significatifs dans ACF a chuté. Nous pouvons enfin commencer la modélisation!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Cours intensif sur les méthodes ARIMA\n",
    "\n",
    "Nous expliquerons ce modèle en construisant lettre par lettre. $SARIMA(p, d, q)(P, D, Q, s)$, (Seasonal Autoregression Moving Average model) modèle de moyenne mobile d'autorégression saisonnière:\n",
    "\n",
    "- $AR(p)$ - modèle d'autorégression, c'est-à-dire régression de la série chronologique sur elle-même. L'hypothèse de base est que les valeurs actuelles de la série dépendent de ses valeurs précédentes avec un certain décalage (ou plusieurs décalages). Le décalage maximum dans le modèle est appelé $p$. Pour déterminer le $p$ initial, vous devez regarder le tracé PACF et trouver le plus grand décalage significatif après lequel **la plupart** des autres décalages deviennent insignifiants.\n",
    "- $MA(q)$ - modèle de moyenne mobile. Sans entrer dans trop de détails, cela modélise l'erreur de la série chronologique, là encore avec l'hypothèse que l'erreur actuelle dépend de la précédente avec un certain retard, ce qui est appelé $q$. La valeur initiale peut être trouvée sur le tracé ACF avec la même logique que précédemment.\n",
    "\n",
    "Combinons nos 4 premières lettres:\n",
    "\n",
    "$AR(p) + MA(q) = ARMA(p, q)$\n",
    "\n",
    "Ce que nous avons ici, c'est le modèle autorégressif à moyenne mobile! Si la série est stationnaire, elle peut être approximée avec ces 4 lettres. Nous allons continuer.\n",
    "\n",
    "- $I(d)$ - ordre d'intégration. Il s'agit simplement du nombre de différences non saisonnières nécessaires pour rendre la série stationnaire. Dans notre cas, c'est juste 1 car nous avons utilisé les premières différences.\n",
    "\n",
    "L'ajout de cette lettre aux quatre autres nous donne le modèle $ARIMA$ qui peut gérer des données non stationnaires à l'aide de différences non saisonnières. Super, une lettre de plus à faire!\n",
    "\n",
    "- $S(s)$ - est responsable de la saisonnalité et est égal à la durée de la saison de la série\n",
    "\n",
    "Avec cela, nous avons trois paramètres: $(P, D, Q)$\n",
    "\n",
    "- $P$ - ordre d'autorégression pour la composante saisonnière du modèle, qui peut être dérivé du PACF. Mais vous devez regarder le nombre de retards importants, qui sont les multiples de la durée de la saison. Par exemple, si la période est égale à 24 et que les décalages des 24e et 48e sont significatifs dans le PACF, cela signifie que le $P$ initial devrait être 2.\n",
    "\n",
    "- $Q$ - logique similaire utilisant le graphique ACF à la place.\n",
    "\n",
    "- $D$ - ordre d'intégration saisonnière. Peut être égal à 1 ou 0, selon que des différences saisonnières ont été appliquées ou non."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Maintenant que nous savons comment définir les paramètres initiaux, regardons à nouveau le tracé final et définissons les paramètres:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsplot(ads_diff[24 + 1 :], lags=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "- $p$ est très probablement 4 puisqu'il s'agit du dernier retard significatif sur le PACF, après quoi la plupart des autres ne sont pas significatifs.\n",
    "- $d$ est égal à 1 car nous avons eu les premières différences\n",
    "- $q$ devrait se situer aux alentours de 4 ainsi que sur l'ACF\n",
    "- $P$ pourrait être égal à 2, car les décalages des 24e et 48e sont quelque peu significatifs sur le PACF\n",
    "- $D$ est à nouveau égal à 1 car nous avons effectué une différenciation saisonnière\n",
    "- $Q$ vaut probablement 1. Le 24 e décalage sur ACF est significatif tandis que le 48 e ne l'est pas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Testons différents modèles et voyons lequel est le meilleur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setting initial values and some bounds for them\n",
    "ps = range(2, 5)\n",
    "d = 1\n",
    "qs = range(2, 5)\n",
    "Ps = range(0, 2)\n",
    "D = 1\n",
    "Qs = range(0, 2)\n",
    "s = 24  # season length is still 24\n",
    "\n",
    "# creating list with all the possible combinations of parameters\n",
    "parameters = product(ps, qs, Ps, Qs)\n",
    "parameters_list = list(parameters)\n",
    "len(parameters_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimizeSARIMA(parameters_list, d, D, s):\n",
    "    \"\"\"\n",
    "        Return dataframe with parameters and corresponding AIC\n",
    "        \n",
    "        parameters_list - list with (p, q, P, Q) tuples\n",
    "        d - integration order in ARIMA model\n",
    "        D - seasonal integration order \n",
    "        s - length of season\n",
    "    \"\"\"\n",
    "\n",
    "    results = []\n",
    "    best_aic = float(\"inf\")\n",
    "\n",
    "    for param in tqdm_notebook(parameters_list):\n",
    "        # we need try-except because on some combinations model fails to converge\n",
    "        try:\n",
    "            model = sm.tsa.statespace.SARIMAX(\n",
    "                ads.Ads,\n",
    "                order=(param[0], d, param[1]),\n",
    "                seasonal_order=(param[2], D, param[3], s),\n",
    "            ).fit(disp=-1)\n",
    "        except:\n",
    "            continue\n",
    "        aic = model.aic\n",
    "        # saving best model, AIC and parameters\n",
    "        if aic < best_aic:\n",
    "            best_model = model\n",
    "            best_aic = aic\n",
    "            best_param = param\n",
    "        results.append([param, model.aic])\n",
    "\n",
    "    result_table = pd.DataFrame(results)\n",
    "    result_table.columns = [\"parameters\", \"aic\"]\n",
    "    # sorting in ascending order, the lower AIC is - the better\n",
    "    result_table = result_table.sort_values(by=\"aic\", ascending=True).reset_index(\n",
    "        drop=True\n",
    "    )\n",
    "\n",
    "    return result_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "result_table = optimizeSARIMA(parameters_list, d, D, s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_table.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the parameters that give the lowest AIC\n",
    "p, q, P, Q = result_table.parameters[0]\n",
    "\n",
    "best_model = sm.tsa.statespace.SARIMAX(\n",
    "    ads.Ads, order=(p, d, q), seasonal_order=(P, D, Q, s)\n",
    ").fit(disp=-1)\n",
    "print(best_model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's inspect the residuals of the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tsplot(best_model.resid[24 + 1 :], lags=60)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Il est clair que les résidus sont stationnaires et il n'y a aucune autocorrélation apparente. Faisons des prédictions en utilisant notre modèle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "def plotSARIMA(series, model, n_steps):\n",
    "    \"\"\"\n",
    "        Plots model vs predicted values\n",
    "        \n",
    "        series - dataset with timeseries\n",
    "        model - fitted SARIMA model\n",
    "        n_steps - number of steps to predict in the future\n",
    "        \n",
    "    \"\"\"\n",
    "    # adding model values\n",
    "    data = series.copy()\n",
    "    data.columns = [\"actual\"]\n",
    "    data[\"arima_model\"] = model.fittedvalues\n",
    "    # making a shift on s+d steps, because these values were unobserved by the model\n",
    "    # due to the differentiating\n",
    "    data[\"arima_model\"][: s + d] = np.NaN\n",
    "\n",
    "    # forecasting on n_steps forward\n",
    "    forecast = model.predict(start=data.shape[0], end=data.shape[0] + n_steps)\n",
    "    forecast = data.arima_model.append(forecast)\n",
    "    # calculate error, again having shifted on s+d steps from the beginning\n",
    "    error = mean_absolute_percentage_error(\n",
    "        data[\"actual\"][s + d :], data[\"arima_model\"][s + d :]\n",
    "    )\n",
    "\n",
    "    plt.figure(figsize=(15, 7))\n",
    "    plt.title(\"Mean Absolute Percentage Error: {0:.2f}%\".format(error))\n",
    "    plt.plot(forecast, color=\"r\", label=\"model\")\n",
    "    plt.axvspan(data.index[-1], forecast.index[-1], alpha=0.5, color=\"lightgrey\")\n",
    "    plt.plot(data.actual, label=\"actual\")\n",
    "    plt.legend()\n",
    "    plt.grid(True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotSARIMA(ads, best_model, 50)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Au final, nous avons obtenu des prévisions très adéquates. Notre modèle s'est trompé de 4,01% en moyenne, ce qui est très, très bon. Cependant, les coûts globaux de préparation des données, de mise à l'arrêt de la série et de sélection des paramètres peuvent ne pas valoir cette précision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Modèles linéaires (et pas tout à fait) pour les séries chronologiques\n",
    "\n",
    "Souvent, dans mon travail, je dois construire des modèles [*rapide, bon, bon marché*](http://fastgood.cheap) comme seul principe directeur. Cela signifie que certains de ces modèles ne seront jamais considérés comme «prêts pour la production» car ils demandent trop de temps pour la préparation des données (comme dans SARIMA) ou nécessitent un recyclage fréquent sur les nouvelles données (encore une fois, SARIMA) ou sont difficiles à régler (bon exemple - SARIMA). Par conséquent, il est très souvent beaucoup plus facile de sélectionner quelques entités à partir des séries chronologiques existantes et de créer un modèle de régression linéaire simple ou, par exemple, une forêt aléatoire. C'est bon et pas cher.\n",
    "\n",
    "Cette approche n'est pas soutenue par la théorie et casse plusieurs hypothèses (par exemple le théorème de Gauss-Markov, en particulier pour les erreurs non corrélées), mais elle est très utile dans la pratique et est souvent utilisée dans les compétitions d'apprentissage automatique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Extraction de caractéristiques"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Le modèle a besoin de caractéristiques, et tout ce que nous avons est une série chronologique à une dimension. Quelles caractéristiques pouvons-nous extraire?\n",
    "\n",
    "* Time series lags (Les décalages des séries chronologiques)\n",
    "* Statistiques de fenêtre:\n",
    "    - Valeur max / min de la série dans une fenêtre\n",
    "    - Valeur moyenne / médiane dans une fenêtre\n",
    "    - Variance de fenêtre\n",
    "    - etc.\n",
    "*Caractéristiques de date et d'heure:\n",
    "    - Minute d'une heure, heure d'un jour, jour de la semaine, etc.\n",
    "    - Ce jour est-il un jour férié? Peut-être qu'il y a un événement spécial? Représentez cela comme une fonction booléenne\n",
    "* Encodage cible\n",
    "* Prévisions d'autres modèles (notez que nous pouvons perdre la vitesse de prédiction de cette façon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Passons en revue certaines des méthodes et voyons ce que nous pouvons extraire de nos données de séries chronologiques sur les annonces.\n",
    "\n",
    "## Retards dans les séries chronologiques (Time series lags)\n",
    "\n",
    "En décalant les pas de la série $n$, nous obtenons une colonne d'entités où la valeur actuelle de la série chronologique est alignée avec sa valeur au moment $t-n$. Si nous effectuons un décalage de 1 et formons un modèle sur cette caractéristique, le modèle sera en mesure de prévoir une longueur d'avance sur l'observation de l'état actuel de la série. Augmenter le décalage, disons jusqu'à 6, permettra au modèle de faire des prévisions 6 pas en avant; cependant, il utilisera les données observées 6 pas en arrière. Si quelque chose change fondamentalement la série au cours de cette période non observée, le modèle n'acceptera pas ces changements et renverra des prévisions avec une grande erreur. Par conséquent, lors de la sélection initiale du décalage, il faut trouver un équilibre entre la qualité de prédiction optimale et la longueur de l'horizon de prévision."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Creating a copy of the initial datagrame to make various transformations\n",
    "data = pd.DataFrame(ads.Ads.copy())\n",
    "data.columns = [\"y\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Adding the lag of the target variable from 6 steps back up to 24\n",
    "for i in range(6, 25):\n",
    "    data[\"lag_{}\".format(i)] = data.y.shift(i)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# take a look at the new dataframe\n",
    "data.tail(7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Très bien, nous avons généré un ensemble de données ici. Pourquoi ne formons-nous pas maintenant un modèle?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "# for time-series cross-validation set 5 folds\n",
    "tscv = TimeSeriesSplit(n_splits=5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def timeseries_train_test_split(X, y, test_size):\n",
    "    \"\"\"\n",
    "        Perform train-test split with respect to time series structure\n",
    "    \"\"\"\n",
    "\n",
    "    # get the index after which test set starts\n",
    "    test_index = int(len(X) * (1 - test_size))\n",
    "\n",
    "    X_train = X.iloc[:test_index]\n",
    "    y_train = y.iloc[:test_index]\n",
    "    X_test = X.iloc[test_index:]\n",
    "    y_test = y.iloc[test_index:]\n",
    "\n",
    "    return X_train, X_test, y_train, y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = data.dropna().y\n",
    "X = data.dropna().drop([\"y\"], axis=1)\n",
    "\n",
    "# reserve 30% of data for testing\n",
    "X_train, X_test, y_train, y_test = timeseries_train_test_split(X, y, test_size=0.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# machine learning in two lines\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0,
     38
    ]
   },
   "outputs": [],
   "source": [
    "def plotModelResults(\n",
    "    model, X_train=X_train, X_test=X_test, plot_intervals=False, plot_anomalies=False\n",
    "):\n",
    "    \"\"\"\n",
    "        Plots modelled vs fact values, prediction intervals and anomalies\n",
    "    \n",
    "    \"\"\"\n",
    "\n",
    "    prediction = model.predict(X_test)\n",
    "\n",
    "    plt.figure(figsize=(15, 7))\n",
    "    plt.plot(prediction, \"g\", label=\"prediction\", linewidth=2.0)\n",
    "    plt.plot(y_test.values, label=\"actual\", linewidth=2.0)\n",
    "\n",
    "    if plot_intervals:\n",
    "        cv = cross_val_score(\n",
    "            model, X_train, y_train, cv=tscv, scoring=\"neg_mean_absolute_error\"\n",
    "        )\n",
    "        mae = cv.mean() * (-1)\n",
    "        deviation = cv.std()\n",
    "\n",
    "        scale = 1.96\n",
    "        lower = prediction - (mae + scale * deviation)\n",
    "        upper = prediction + (mae + scale * deviation)\n",
    "\n",
    "        plt.plot(lower, \"r--\", label=\"upper bond / lower bond\", alpha=0.5)\n",
    "        plt.plot(upper, \"r--\", alpha=0.5)\n",
    "\n",
    "        if plot_anomalies:\n",
    "            anomalies = np.array([np.NaN] * len(y_test))\n",
    "            anomalies[y_test < lower] = y_test[y_test < lower]\n",
    "            anomalies[y_test > upper] = y_test[y_test > upper]\n",
    "            plt.plot(anomalies, \"o\", markersize=10, label=\"Anomalies\")\n",
    "\n",
    "    error = mean_absolute_percentage_error(prediction, y_test)\n",
    "    plt.title(\"Mean absolute percentage error {0:.2f}%\".format(error))\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.tight_layout()\n",
    "    plt.grid(True)\n",
    "\n",
    "\n",
    "def plotCoefficients(model):\n",
    "    \"\"\"\n",
    "        Plots sorted coefficient values of the model\n",
    "    \"\"\"\n",
    "\n",
    "    coefs = pd.DataFrame(model.coef_, X_train.columns)\n",
    "    coefs.columns = [\"coef\"]\n",
    "    coefs[\"abs\"] = coefs.coef.apply(np.abs)\n",
    "    coefs = coefs.sort_values(by=\"abs\", ascending=False).drop([\"abs\"], axis=1)\n",
    "\n",
    "    plt.figure(figsize=(15, 7))\n",
    "    coefs.coef.plot(kind=\"bar\")\n",
    "    plt.grid(True, axis=\"y\")\n",
    "    plt.hlines(y=0, xmin=0, xmax=len(coefs), linestyles=\"dashed\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotModelResults(lr, plot_intervals=True)\n",
    "plotCoefficients(lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Des retards simples et une régression linéaire nous ont donné des prédictions qui ne sont pas si éloignées de SARIMA en termes de qualité. Il existe de nombreuses caractéristiques inutiles, nous ferons donc la sélection des caractéristiques dans un petit moment. Pour l'instant, continuons l'ingénierie!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous ajouterons l'heure, le jour de la semaine et un booléen pour `is_weekend`. Pour ce faire, nous devons transformer l'index actuel de la trame de données au format `datetime` et extraire` hour` et `week`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.index = pd.to_datetime(data.index)\n",
    "data[\"hour\"] = data.index.hour\n",
    "data[\"weekday\"] = data.index.weekday\n",
    "data[\"is_weekend\"] = data.weekday.isin([5, 6]) * 1\n",
    "data.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous pouvons visualiser les caractéristiques résultantes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16, 5))\n",
    "plt.title(\"Encoded features\")\n",
    "data.hour.plot()\n",
    "data.weekday.plot()\n",
    "data.is_weekend.plot()\n",
    "plt.grid(True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Étant donné que nous avons maintenant différentes échelles dans nos variables, des milliers pour les caractéristiques de décalage et des dizaines pour les caractéristiques catégorielles, nous devons les transformer en la même échelle pour explorer l'importance des caractéristiques et, plus tard, la régularisation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "\n",
    "scaler = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = data.dropna().y\n",
    "X = data.dropna().drop([\"y\"], axis=1)\n",
    "\n",
    "X_train, X_test, y_train, y_test = timeseries_train_test_split(X, y, test_size=0.3)\n",
    "\n",
    "X_train_scaled = scaler.fit_transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)\n",
    "\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train_scaled, y_train)\n",
    "\n",
    "plotModelResults(lr, X_train=X_train_scaled, X_test=X_test_scaled, plot_intervals=True)\n",
    "plotCoefficients(lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "L'erreur de test diminue un peu. A en juger par le tracé des coefficients, nous pouvons dire que `week` et` is_weekend` sont des caractéristiques utiles."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Encodage de la  cible\n",
    "Je voudrais ajouter une autre variante pour l'encodage des variables catégorielles: l'encodage par valeur moyenne. S'il n'est pas souhaitable d'exploser un ensemble de données en utilisant de nombreuses variables factices qui peuvent entraîner la perte d'informations et si elles ne peuvent pas être utilisées comme valeurs réelles en raison de conflits comme \"0 heures <23 heures\", il est alors possible de coder une variable avec des valeurs légèrement plus interprétables. L'idée naturelle est de coder avec la valeur moyenne de la variable cible. Dans notre exemple, chaque jour de la semaine et chaque heure de la journée peut être codé par le nombre moyen correspondant d'annonces vues au cours de cette journée ou de cette heure. Il est très important de s'assurer que la valeur moyenne est calculée sur l'ensemble d'apprentissage uniquement (ou sur le 'pli' (fold) de validation croisée actuel uniquement) afin que le modèle ne soit pas au courant de l'avenir."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def code_mean(data, cat_feature, real_feature):\n",
    "    \"\"\"\n",
    "    Returns a dictionary where keys are unique categories of the cat_feature,\n",
    "    and values are means over real_feature\n",
    "    \"\"\"\n",
    "    return dict(data.groupby(cat_feature)[real_feature].mean())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Regardons les moyennes par heure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "average_hour = code_mean(data, \"hour\", \"y\")\n",
    "plt.figure(figsize=(7, 5))\n",
    "plt.title(\"Hour averages\")\n",
    "pd.DataFrame.from_dict(average_hour, orient=\"index\")[0].plot()\n",
    "plt.grid(True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Enfin, rassemblons toutes les transformations dans une seule fonction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "def prepareData(series, lag_start, lag_end, test_size, target_encoding=False):\n",
    "    \"\"\"\n",
    "        series: pd.DataFrame\n",
    "            dataframe with timeseries\n",
    "\n",
    "        lag_start: int\n",
    "            initial step back in time to slice target variable \n",
    "            example - lag_start = 1 means that the model \n",
    "                      will see yesterday's values to predict today\n",
    "\n",
    "        lag_end: int\n",
    "            final step back in time to slice target variable\n",
    "            example - lag_end = 4 means that the model \n",
    "                      will see up to 4 days back in time to predict today\n",
    "\n",
    "        test_size: float\n",
    "            size of the test dataset after train/test split as percentage of dataset\n",
    "\n",
    "        target_encoding: boolean\n",
    "            if True - add target averages to the dataset\n",
    "        \n",
    "    \"\"\"\n",
    "\n",
    "    # copy of the initial dataset\n",
    "    data = pd.DataFrame(series.copy())\n",
    "    data.columns = [\"y\"]\n",
    "\n",
    "    # lags of series\n",
    "    for i in range(lag_start, lag_end):\n",
    "        data[\"lag_{}\".format(i)] = data.y.shift(i)\n",
    "\n",
    "    # datetime features\n",
    "    data.index = pd.to_datetime(data.index)\n",
    "    data[\"hour\"] = data.index.hour\n",
    "    data[\"weekday\"] = data.index.weekday\n",
    "    data[\"is_weekend\"] = data.weekday.isin([5, 6]) * 1\n",
    "\n",
    "    if target_encoding:\n",
    "        # calculate averages on train set only\n",
    "        test_index = int(len(data.dropna()) * (1 - test_size))\n",
    "        data[\"weekday_average\"] = list(\n",
    "            map(code_mean(data[:test_index], \"weekday\", \"y\").get, data.weekday)\n",
    "        )\n",
    "        data[\"hour_average\"] = list(\n",
    "            map(code_mean(data[:test_index], \"hour\", \"y\").get, data.hour)\n",
    "        )\n",
    "\n",
    "        # drop encoded variables\n",
    "        data.drop([\"hour\", \"weekday\"], axis=1, inplace=True)\n",
    "\n",
    "    # train-test split\n",
    "    y = data.dropna().y\n",
    "    X = data.dropna().drop([\"y\"], axis=1)\n",
    "    X_train, X_test, y_train, y_test = timeseries_train_test_split(\n",
    "        X, y, test_size=test_size\n",
    "    )\n",
    "\n",
    "    return X_train, X_test, y_train, y_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = prepareData(\n",
    "    ads.Ads, lag_start=6, lag_end=25, test_size=0.3, target_encoding=True\n",
    ")\n",
    "\n",
    "X_train_scaled = scaler.fit_transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)\n",
    "\n",
    "lr = LinearRegression()\n",
    "lr.fit(X_train_scaled, y_train)\n",
    "\n",
    "plotModelResults(\n",
    "    lr,\n",
    "    X_train=X_train_scaled,\n",
    "    X_test=X_test_scaled,\n",
    "    plot_intervals=True,\n",
    "    plot_anomalies=True,\n",
    ")\n",
    "plotCoefficients(lr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous voyons un **surajustement**! `Hour_average` était si bien dans l'ensemble de données d'entraînement que le modèle a décidé de concentrer toutes ses forces sur celui-ci. En conséquence, la qualité de la prévision a chuté. Ce problème peut être résolu de différentes manières; par exemple, nous pouvons calculer l'encodage cible non pas pour l'ensemble du train, mais pour une certaine fenêtre à la place. De cette façon, les encodages de la dernière fenêtre observée décrivent probablement mieux l'état actuel de la série. Alternativement, nous pouvons simplement le supprimer manuellement car nous sommes sûrs que cela ne fait qu'empirer les choses dans ce cas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = prepareData(\n",
    "    ads.Ads, lag_start=6, lag_end=25, test_size=0.3, target_encoding=False\n",
    ")\n",
    "\n",
    "X_train_scaled = scaler.fit_transform(X_train)\n",
    "X_test_scaled = scaler.transform(X_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "## Régularisation et sélection des caractéristiques\n",
    "\n",
    "Comme nous le savons déjà, toutes les fonctionnalités ne sont pas également saines - certaines peuvent entraîner un sur-ajustement tandis que d'autres doivent être supprimées. Outre l'inspection manuelle, nous pouvons appliquer la régularisation. Les régressions Ridge et Lasso sont deux des modèles de régression les plus populaires avec régularisation. Ils ajoutent tous deux des contraintes supplémentaires à notre fonction de perte.\n",
    "\n",
    "Dans le cas de la régression Ridge, ces contraintes sont la somme des carrés des coefficients multipliés par le coefficient de régularisation. Plus le coefficient d'une caractéristique est élevé, plus notre perte sera importante. Par conséquent, nous essaierons d'optimiser le modèle tout en maintenant les coefficients assez bas.\n",
    "\n",
    "À la suite de cette régularisation $L2$, nous aurons un biais plus élevé et une variance plus faible, donc le modèle généralisera mieux (du moins c'est ce que nous espérons arriver).\n",
    "\n",
    "Le deuxième modèle de régression, la régression Lasso, ajoute à la fonction de perte, non des carrés, mais des valeurs absolues des coefficients. Par conséquent, au cours du processus d'optimisation, les coefficients des entités non importantes peuvent devenir des zéros, ce qui permet une sélection automatisée des entités. Ce type de régularisation est appelé $L1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Tout d'abord, assurons-nous que nous avons des caractéristiques à supprimer et que les données ont des caractéristiques hautement corrélées."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 8))\n",
    "sns.heatmap(X_train.corr());"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LassoCV, RidgeCV\n",
    "\n",
    "ridge = RidgeCV(cv=tscv)\n",
    "ridge.fit(X_train_scaled, y_train)\n",
    "\n",
    "plotModelResults(\n",
    "    ridge,\n",
    "    X_train=X_train_scaled,\n",
    "    X_test=X_test_scaled,\n",
    "    plot_intervals=True,\n",
    "    plot_anomalies=True,\n",
    ")\n",
    "plotCoefficients(ridge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous pouvons clairement voir que certains coefficients se rapprochent de plus en plus de zéro (bien qu'ils ne l'atteignent jamais réellement) à mesure que leur importance dans le modèle diminue."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lasso = LassoCV(cv=tscv)\n",
    "lasso.fit(X_train_scaled, y_train)\n",
    "\n",
    "plotModelResults(\n",
    "    lasso,\n",
    "    X_train=X_train_scaled,\n",
    "    X_test=X_test_scaled,\n",
    "    plot_intervals=True,\n",
    "    plot_anomalies=True,\n",
    ")\n",
    "plotCoefficients(lasso)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "La régression Lasso s'est avérée plus conservatrice; il a supprimé le 23 e décalage des caractéristiques les plus importantes et a complètement supprimé 5 caractéristiques, ce qui n'a fait qu'améliorer la qualité des prévisions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "en"
   },
   "source": [
    "# Boosting \n",
    "Why shouldn't we try XGBoost now?\n",
    "<img src=\"../../img/xgboost_the_things.jpg\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Boosting\n",
    "Pourquoi ne devrions-nous pas essayer XGBoost maintenant?\n",
    "<img src=\"https://raw.githubusercontent.com/Yorko/mlcourse.ai/master/img/xgboost-all-the-things.jpg\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     3
    ]
   },
   "outputs": [],
   "source": [
    "from xgboost import XGBRegressor\n",
    "\n",
    "xgb = XGBRegressor(verbosity=0)\n",
    "xgb.fit(X_train_scaled, y_train);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plotModelResults(\n",
    "    xgb,\n",
    "    X_train=X_train_scaled,\n",
    "    X_test=X_test_scaled,\n",
    "    plot_intervals=True,\n",
    "    plot_anomalies=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "Nous avons un gagnant! Il s'agit de la plus petite erreur sur l'ensemble de test parmi tous les modèles que nous avons essayés jusqu'à présent.\n",
    "\n",
    "Mais, cette victoire est trompeuse, et ce n'est peut-être pas l'idée la plus brillante d'intégrer `xgboost` dès que vous obtenez les données de séries chronologiques. En règle générale, les modèles basés sur les arbres gèrent mal les tendances des données par rapport aux modèles linéaires. Dans ce cas, vous devrez d'abord détruire votre série ou utiliser des astuces pour que la magie se produise. Idéalement, vous pouvez rendre la série stationnaire, puis utiliser XGBoost. Par exemple, vous pouvez prévoir la tendance séparément avec un modèle linéaire, puis ajouter des prévisions à partir de `xgboost` pour obtenir une prévision finale."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Conclusion\n",
    "\n",
    "Nous avons discuté de différentes méthodes d'analyse et de prévision des séries chronologiques. Malheureusement, ou peut-être heureusement, il n'y a pas une seule façon de résoudre ce genre de problèmes. Les méthodes développées dans les années 1960 (et certaines même au début du XXIe siècle) sont toujours populaires, tout comme les LSTM et les RNN (non couvertes dans cet article). Cela est partiellement lié au fait que la tâche de prédiction, comme toute autre tâche liée aux données, nécessite de la créativité à bien des égards et nécessite certainement des recherches. Malgré le grand nombre de métriques de qualité formelles et d'approches d'estimation des paramètres, il est souvent nécessaire d'essayer quelque chose de différent pour chaque série chronologique. Enfin et surtout, l'équilibre entre qualité et coût est important. À titre d'exemple, le modèle SARIMA peut produire des résultats spectaculaires après réglage, mais peut nécessiter de nombreuses heures de ~~danse du tambourin~~ manipulation de séries temporelles tandis qu'un modèle de régression linéaire simple peut être construit en 10 minutes et peut obtenir des résultats plus ou moins comparables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lang": "fr"
   },
   "source": [
    "# Ressources utiles\n",
    "\n",
    "* The same notebook as an interactive web-based [Kaggle Kernel](https://www.kaggle.com/kashnitsky/topic-9-part-1-time-series-analysis-in-python)\n",
    "* \"LSTM (Long Short Term Memory) Networks for predicting Time Series\" - a tutorial by Max Sergei Bulaev within mlcourse.ai (full list of tutorials is [here](https://mlcourse.ai/tutorials))\n",
    "* Main course [site](https://mlcourse.ai), [course repo](https://github.com/Yorko/mlcourse.ai), and YouTube [channel](https://www.youtube.com/watch?v=QKTuw4PNOsU&list=PLVlY_7IJCMJeRfZ68eVfEcu-UcN9BbwiX)\n",
    "* Course materials as a [Kaggle Dataset](https://www.kaggle.com/kashnitsky/mlcourse)\n",
    "* Medium [\"story\"](https://medium.com/open-machine-learning-course/open-machine-learning-course-topic-9-time-series-analysis-in-python-a270cb05e0b3?source=collection_home---6------2---------------------) based on this notebook\n",
    "* If you read Russian: an [article](https://habr.com/ru/company/ods/blog/327242/) on Habr.com with ~ the same material. And a [lecture](https://youtu.be/_9lBwXnbOd8) on YouTube\n",
    "* [Online textbook](https://people.duke.edu/~rnau/411home.htm) for the advanced statistical forecasting course at Duke University - covers various smoothing techniques in detail along with linear and ARIMA models\n",
    "* [Comparison of ARIMA and Random Forest time series models for prediction of avian influenza H5N1 outbreaks](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-15-276) - one of a few cases where using random forest for time series forecasting is actively defended\n",
    "* [Time Series Analysis (TSA) in Python - Linear Models to GARCH](http://www.blackarbs.com/blog/time-series-analysis-in-python-linear-models-to-garch/11/1/2016) - applying the ARIMA models family to the task of modeling financial indicators (by Brian Christopher)Christopher)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "nbTranslate": {
   "displayLangs": [
    "*"
   ],
   "hotkey": "alt-t",
   "langInMainMenu": true,
   "sourceLang": "en",
   "targetLang": "fr",
   "useGoogleTranslate": true
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
